library(readxl)
library(httr)
library(arules)
library(tidyr)
library(dplyr)
library(lubridate)
library(stringr)
library(ggplot2)
library(purrr)
library(broom)
library(tibble)
library(plotly)
Odczyt danych przeprowdzony jest bezośrednio ze strony, aby zapewenić stałą wartość początkową dla każdego uruchomienia. Dodatkowo uzupełniane są brakujące dane identyfikatory pacjentów PATIENT_ID. Wartości brakujące wynikają ze sposobu przetwarzania scalonych komórek pliku źródłowego.
# https://stackoverflow.com/a/41368947
url <- "http://www.cs.put.poznan.pl/dbrzezinski/teaching/zed/wuhan_blood_sample_data_Jan_Feb_2020.xlsx"
GET(url, write_disk(tf <- tempfile(fileext = ".xlsx")))
df <- read_excel(tf) %>% fill(PATIENT_ID)
Surowy zbiór danych zawiera 81 atrybutów i 6120 wierszy.
Atrybuty zbioru danych można potwierdzić na dwie grupy:
Poza danymi niedostępnymi (NA) w zbiorze danych wsytępują dane brakujące, oznakowane wartością -1. Poniższa tabela pokazuje ilość takich wartości w każdej kolumnie, o ile te wartości występują.
df %>% pivot_longer(-c(RE_DATE,`Admission time`, `Discharge time`)) %>% filter(value==-1) %>% count(name)
## # A tibble: 2 x 2
## name n
## <chr> <int>
## 1 2019-nCoV nucleic acid detection 501
## 2 Platelet count 4
Dla kolumny Platelet count warkości brakujących jest mało, liczby -1 zostaną zamienione na NA.
df <- df %>% mutate(`Platelet count` = if(is.na(`Platelet count`) || `Platelet count` != -1) `Platelet count` else NA)
Dla kolumny 2019-nCoV nucleic acid detection sprawdzimy całkowitą ilość wartości różnych od NA.
df %>%
pivot_longer(-c(RE_DATE,`Admission time`, `Discharge time`)) %>%
filter(name=="2019-nCoV nucleic acid detection") %>%
filter(!is.na(value)) %>% count(name)
## # A tibble: 1 x 2
## name n
## <chr> <int>
## 1 2019-nCoV nucleic acid detection 501
Wszystkie wartości w kolumnie 2019-nCoV nucleic acid detection różne od NA wynoszą -1, w związku z tym ta kolumna jest właściwie nieprzydatna. Mimo to na ten moment postanowiono nie zmieniać wartości w tej kolumnie, ewentualne zmiany można wykonać w dalszych krokach.
Pierwszych siedem kolumn to atrybuty ogólne, reprezentują następująco:
Nazwy wyżej wymienionych atrybutów zostaną przeformatowane, będą pisane wielkimi literami, a ewentualne spacje będą zamieniane na znak podkreślinka _, w ten sposób atrbuty ogólne będą łatwo odróżnialne od pozostałych atrybutów.
df <- df %>% rename_with(~ toupper(gsub(" ", "_", .x, fixed = TRUE)), 1:7)
Pola GENDER i OUTCOME zostaną dodatkowo zmienione na typ factor - z typu numerycznego zamienione zostaną na zmienne nominalne. Do tego trzeba przeprowadzić identyfikację typu płci.
Wartości płci nie są jawnie opisane, w artykule źródłowym została natomiast podana proporcja pacjentów danej płci, co powinno umożliwić okreslenie właściwego przypisania.
knitr::kable(
df %>% group_by(PATIENT_ID) %>% summarise_all(last) %>% count(GENDER) %>% rename("count" = "n"),
align="lc"
)
| GENDER | count |
|---|---|
| 1 | 224 |
| 2 | 151 |
Wniosek:
12df <- df %>% mutate(across("GENDER", ~factor(., levels=c(1,2), labels=c("MALE", "FEMALE")))) %>%
mutate(across("OUTCOME", ~factor(., levels=c(0,1), labels=c("CURED","DECEASED"))))
Poniżej można zapoznać się z analizą wartości dla atrybutów ogólnych. Wartość atrybutu RE_DATE ma unikalną wartość dla kazdego wiersza.
knitr::kable(summary(select(df, RE_DATE)))
| RE_DATE | |
|---|---|
| Min. :2020-01-10 19:45:00 | |
| 1st Qu.:2020-02-04 13:44:00 | |
| Median :2020-02-09 12:42:30 | |
| Mean :2020-02-08 07:00:02 | |
| 3rd Qu.:2020-02-13 10:34:00 | |
| Max. :2020-02-18 17:49:00 | |
| NA’s :14 |
Pozostałe z atrybtów ogólnych, są wspólne dla każdego pacjenta, zastosowano transformację do zniwelowania skutków duplikacji.
options(knitr.kable.NA = '')
knitr::kable(
summary(
df %>%
select(c(1,3:7)) %>%
group_by(PATIENT_ID) %>%
summarise_all(last)
)
)
| PATIENT_ID | AGE | GENDER | ADMISSION_TIME | DISCHARGE_TIME | OUTCOME | |
|---|---|---|---|---|---|---|
| Min. : 1.0 | Min. :18.00 | MALE :224 | Min. :2020-01-10 15:52:20 | Min. :2020-01-23 09:09:23 | CURED :201 | |
| 1st Qu.: 94.5 | 1st Qu.:46.00 | FEMALE:151 | 1st Qu.:2020-02-01 19:27:40 | 1st Qu.:2020-02-11 13:39:21 | DECEASED:174 | |
| Median :188.0 | Median :62.00 | Median :2020-02-04 22:30:34 | Median :2020-02-16 17:40:07 | |||
| Mean :188.0 | Mean :58.83 | Mean :2020-02-04 20:13:51 | Mean :2020-02-15 16:42:59 | |||
| 3rd Qu.:281.5 | 3rd Qu.:70.00 | 3rd Qu.:2020-02-10 04:11:10 | 3rd Qu.:2020-02-19 11:47:14 | |||
| Max. :375.0 | Max. :95.00 | Max. :2020-02-17 21:30:07 | Max. :2020-03-04 16:21:51 |
Liczba unikatowych identyfikatorów pacjentów wynosi 375.
TODO
Ta sekcja poświęcona jest analizie długości hospitalizacji poszczególnych grup pacjentów.
| HOSPITALIZATION_TIME | quantile |
|---|---|
| 1 days | min |
| 6 days | 1st |
| 11 days | 2nd |
| 17 days | 3rd |
| 36 days | max |
| OUTCOME | GENDER | AGE | HOSPITALIZATION_TIME |
|---|---|---|---|
| CURED | MALE | 54.5 | 15.0 days |
| CURED | FEMALE | 45.0 | 15.0 days |
| DECEASED | MALE | 68.0 | 7.0 days |
| DECEASED | FEMALE | 72.0 | 6.5 days |
Pozostałe 74 kolumn zawiera wartości atrybutów, które reprezentują rodzaje miar zebranych podczas każdego badania.
Nazwy kolumn zostały poddane następującym transformacjom:
(%) zamieniono na suffix _percent(#) usunięto i myślnika - zamieniono na podkreślenie dolne _ ze względu na łatwość odwołania się do kolumnPoniższy blok kodu szczegółowo ukazuje zatosowane transfomrmacje, jednocześnie prezentowana są ostatenczne nazwy kolumn wyników badań.
df <- df %>%
rename("NT-proBNP" = "Amino-terminal brain natriuretic peptide precursor(NT-proBNP)") %>%
rename("2019-nCov-detection" = "2019-nCoV nucleic acid detection") %>%
rename("aPPT" = "Activation of partial thromboplastin time") %>%
rename("ALP" = "Alkaline phosphatase") %>%
rename("AAT" = "aspartate aminotransferase") %>%
rename("FDPs" = "Fibrin degradation products") %>%
rename("ALT" = "glutamic-pyruvic transaminase") %>%
rename("hs-CRP" = "High sensitivity C-reactive protein") %>%
rename("hs-cTnI" = "Hypersensitive cardiac troponinI") %>%
rename("HCV-abs-quant" = "HCV antibody quantification") %>%
rename("HIV-abs-quant" = "HIV antibody quantification") %>%
rename("INR" = "International standard ratio") %>%
rename("LDH" = "Lactate dehydrogenase") %>%
rename("MCH" = "mean corpuscular hemoglobin") %>%
rename("MCHC" = "mean corpuscular hemoglobin concentration") %>%
rename("MCV" = "mean corpuscular volume") %>%
rename("MPV" = "Mean platelet volume") %>%
rename("P-LCR" = "platelet large cell ratio") %>%
rename("PDW" = "PLT distribution width") %>%
rename("q-t-pallidum-abs" = "Quantification of Treponema pallidum antibodies") %>%
rename("RCDW-SD" = "RBC distribution width SD") %>%
rename("RBC_count" = "Red blood cell count") %>%
rename("RCDW" = "Red blood cell distribution width") %>%
rename("TNF-alfa" = "Tumor necrosis factorα") %>%
rename("WBC_count" = "White blood cell count") %>%
rename("gamma-GT" = "γ-glutamyl transpeptidase") %>%
rename_with(~str_c(str_replace(.x, "\\(%\\)", ""), "_percent"), contains("(%)")) %>%
rename_with(~str_replace(.x, "\\(#\\)", "")) %>%
rename_with(tolower, matches("^[[:upper:]]{1}[[:lower:]]+", ignore.case = FALSE)) %>%
rename_with(~str_replace_all(.x, "[ |-]", "_"))
str_sort(colnames(df[,-(1:7)]))
## [1] "2019_nCov_detection" "AAT" "albumin"
## [4] "ALP" "ALT" "antithrombin"
## [7] "aPPT" "basophil_count" "basophil_percent"
## [10] "calcium" "corrected_calcium" "creatinine"
## [13] "D_D_dimer" "direct_bilirubin" "eGFR"
## [16] "eosinophil_count" "eosinophils_percent" "ESR"
## [19] "FDPs" "ferritin" "fibrinogen"
## [22] "gamma_GT" "globulin" "glucose"
## [25] "HBsAg" "HCO3_" "HCV_abs_quant"
## [28] "hematocrit" "hemoglobin" "HIV_abs_quant"
## [31] "hs_CRP" "hs_cTnI" "indirect_bilirubin"
## [34] "INR" "interleukin_10" "interleukin_1β"
## [37] "interleukin_2_receptor" "interleukin_6" "interleukin_8"
## [40] "LDH" "lymphocyte_count" "lymphocyte_percent"
## [43] "MCH" "MCHC" "MCV"
## [46] "monocytes_count" "monocytes_percent" "MPV"
## [49] "neutrophils_count" "neutrophils_percent" "NT_proBNP"
## [52] "P_LCR" "PDW" "PH_value"
## [55] "platelet_count" "procalcitonin" "prothrombin_activity"
## [58] "prothrombin_time" "q_t_pallidum_abs" "RBC_count"
## [61] "RCDW" "RCDW_SD" "serum_chloride"
## [64] "serum_potassium" "serum_sodium" "thrombin_time"
## [67] "thrombocytocrit" "TNF_alfa" "total_bilirubin"
## [70] "total_cholesterol" "total_protein" "urea"
## [73] "uric_acid" "WBC_count"
Niniejsza sekcja jest poświęcona wstępnemu przetwarzniu zbiorów danych, ze szczególnym uwzględnieniem usuwania i scalania wybranych rekordów.
W powyższej głównej zauważono wystąpienia wartości brakujących dla kolumny RE_DATE. Szybkie sprawdzenie wykazuje, że dla wierszy, gdzie nie ma atrybut RE_DATE przyjmuje wartość ‘NA’, wartości wszystkich atrybutów pochodzących z badań krwi, również przyjmują wartość ‘NA’.
all(sapply(df[is.na(df$RE_DATE), -c(1,3:7)], is.na))
## [1] TRUE
Fakt, że wiersze te są właściwie puste, jest podstawą do ich usunięcia ze zbioru.
df <- df[!is.na(df$RE_DATE),]
Po tej zmianie w zbiorze pozostało 6106 wierszy i 361 unikalnych identyfikatorów pacjentów.
Wstępnę analiza danych sugeruje wskazuje, że wiele wartości atrybutów jest oznaczona wartością ‘NA’. Poniższy wykresy wskazuje, ilość wystąpień rekordów badań z danę liczbą wartości (różnych od NA).
Z powyższego histogramu można wyciągnąć kilka wniosków
W poprzedniej sekcji stwierdzono, że można przeanalizować współwystępowanie ze sobą poszczególnych atrybutów w rekordach. W związku z tym tabela wyników badań zostanie przekształcona do zbioru binarnych transakcji.
transactions <- sapply(df[,-(1:7)], function(x) !is.na(x))
Sprawdzona zostanie ilość unikalnych transakcji.
length(unique(apply(transactions, 1, function(x) which(x))))
## [1] 176
W zbiorze 6106 transakcji jest zaledwie 176 różnych typów.
Aby sprawdzić, czy atrybuty wsytępują w grupach, wyszukiwane są maksymalne domknięte zbiory częste, dla mimalnej wartości support równej 5%.
itemsets <- apriori(transactions, parameter = list(target="closed frequent itemsets", support=.05, minlen=1, maxlen=30, maxtime=60))
## Apriori
##
## Parameter specification:
## confidence minval smax arem aval originalSupport maxtime support minlen
## NA 0.1 1 none FALSE TRUE 60 0.05 1
## maxlen target ext
## 30 closed frequent itemsets TRUE
##
## Algorithmic control:
## filter tree heap memopt load sort verbose
## 0.1 TRUE TRUE FALSE TRUE 2 TRUE
##
## Absolute minimum support count: 305
##
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[74 item(s), 6106 transaction(s)] done [0.01s].
## sorting and recoding items ... [63 item(s)] done [0.00s].
## creating transaction tree ... done [0.00s].
## checking subsets of size 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 done [27.43s].
## filtering closed item sets ... done [15.95s].
## sorting transactions ... done [0.00s].
## writing ... [130 set(s)] done [0.34s].
## creating S4 object ... done [1.02s].
mc_itemsets <- itemsets[is.maximal(itemsets)]
inspect(sort(mc_itemsets, by = 'support'))
## items support transIdenticalToItemsets count
## [1] {hemoglobin,
## eosinophils_percent,
## basophil_percent,
## platelet_count,
## monocytes_percent,
## RCDW,
## neutrophils_percent,
## MCV,
## hematocrit,
## WBC_count,
## MCHC,
## lymphocyte_count,
## RBC_count,
## eosinophil_count,
## neutrophils_count,
## MPV,
## RCDW_SD,
## lymphocyte_percent,
## P_LCR,
## monocytes_count,
## PDW,
## basophil_count,
## MCH,
## thrombocytocrit} 0.13609564 0.1354406 831
## [2] {glucose} 0.12692434 0.0000000 775
## [3] {serum_chloride,
## ALP,
## albumin,
## total_bilirubin,
## indirect_bilirubin,
## total_protein,
## urea,
## corrected_calcium,
## serum_potassium,
## direct_bilirubin,
## total_cholesterol,
## AAT,
## uric_acid,
## HCO3_,
## calcium,
## LDH,
## globulin,
## gamma_GT,
## hs_CRP,
## serum_sodium,
## ALT,
## eGFR,
## creatinine} 0.11021946 0.1094006 673
## [4] {hs_cTnI} 0.08303308 0.0000000 507
## [5] {2019_nCov_detection} 0.08205044 0.0000000 501
## [6] {NT_proBNP} 0.07779234 0.0000000 475
## [7] {procalcitonin} 0.07517196 0.0000000 459
## [8] {PH_value} 0.06288896 0.0000000 384
## [9] {ESR} 0.06272519 0.0000000 383
## [10] {prothrombin_time,
## antithrombin,
## prothrombin_activity,
## fibrinogen,
## thrombin_time,
## D_D_dimer,
## FDPs,
## INR,
## aPPT} 0.05371765 0.0000000 328
Otrzymaliśmy 10 maksymalnych domkniętych zbiorów częstych, z czego 7 to zbiory jednolementowe. Pozostałe 3 zbiory (1, 3 i 10) zawierają wiele elementów. Na podstawie elementów składowych, zbiory wieolemelemnetowe zostały zakwalifikowane następująco:
Na podstawie uzyskanego wyniku można wnioskować, że ilość elementów w rekordzie może wynikać z procesu technologicznego analizowania próbek krwi. Sugeruje to, że możliwe będzie scalanie rekordów, np na podstawie dnia wykonania badania (jeśli jednemu pacjentowi przeprowadzono wiele badań w ciagu dnia).
W tej sekcji przeprowadzona zostanie analiza występowania badań w czasie. Ponadto rekordy zostaną pogrupowane według przynależności do zbiorów wyznaczonych w poprzedniej sekcji oraz dwóch dodatkowych klas NONE jeśli rekord nie należy do żadnej kategorii i MULTIPLE jeśli rekod należy do więcej niż jednej kategorii.
Poniższe blok kodu generuje oetykietowany zbiór rekordów badań, ograniczony jedynie do podstawowych informacji.
class_labels <- as(items(sort(mc_itemsets, by="support")), "list")
class_labels <- sapply(class_labels, function(x) x[1])
class_labels[1] <- "Complete blood count"
class_labels[3] <- "Blood biochemistry"
class_labels[10] <- "Coagulation"
classes <- as(items(sort(mc_itemsets, by="support")), "matrix")
labeled <- df %>% rowwise() %>%
mutate(CLASS = paste(
which(
apply(
(matrix(
!is.na(c_across(-(1:7))),
nrow=nrow(classes),
ncol=ncol(classes),
byrow=TRUE,
dimnames=dimnames(classes)
) & classes) == classes,
1,
all)
), collapse=""),
.after=OUTCOME ) %>%
mutate(CLASS = if(CLASS != "" && as.integer(CLASS)>10) "MULTIPLE" else CLASS) %>%
mutate(CLASS = if(CLASS != "" && CLASS != "MULTIPLE") class_labels[as.integer(CLASS)] else CLASS) %>%
mutate(CLASS = str_replace(CLASS, "^$", "NONE")) %>%
mutate(CLASS = str_trim(CLASS)) %>%
mutate(CLASS = factor(CLASS)) %>%
rowwise() %>%
mutate(SIZE = sum(!is.na(c_across(-c(1:8)))), .after=CLASS) %>%
select(1:9) %>%
ungroup()
knitr::kable(
labeled %>% group_by(CLASS) %>% summarize(count = n(), .groups="drop") %>% arrange(desc(count))
)
| CLASS | count |
|---|---|
| NONE | 1433 |
| Complete blood count | 816 |
| glucose | 584 |
| Blood biochemistry | 507 |
| 2019_nCov_detection | 500 |
| MULTIPLE | 485 |
| hs_cTnI | 386 |
| PH_value | 375 |
| ESR | 368 |
| Coagulation | 300 |
| NT_proBNP | 190 |
| procalcitonin | 162 |
Powyższa tabela prezentuje ilość rekordów przydzielonych do każdej z klas.
Poniższy wykres prezentuje badania danej kategorii według daty wykonania. Dodatkowo każdy punkt jest opisany identyfikatorem pacjenta.
Z powyższego wykresu można wywnioskować kilka rzeczy:
RE_DATE. Biorąc pod uwagę fakt, że takie punkty najcześciej przynależą do różnych pacjentów i do tej samej kategorii, należy wnioskować, że poszczególne rodzaje badań były wykonywane równolegle dla wielu pacjentów.Jednocześnie nie znaleziono żadnych czynników, które podważają sens scalania rekordów.
TODO Analiza ilości badań przeprowadzonych jednemu pacjentowi.
patient_last_results <- df %>%
group_by(PATIENT_ID) %>%
fill(everything()) %>%
summarise_all(last) %>%
rowwise() %>%
mutate(number_of_tests=sum(!is.na(c_across(-(1:7)))), .after=OUTCOME)
ggplot(
patient_last_results,
aes(number_of_tests)
) + geom_histogram(binwidth = 1)
Zauważamy, że część pacjentów ma niewiele przeprowadzonych rodzajów badań.
TODO analityczne wybranie tej wartości
veryfication_df <- function(n_tests, cols_cnt, th_seq) {
cls_vec <- c()
pts_vec <- c()
for(th in th_seq) {
tmp <- n_tests %>% filter(number_of_tests < cols_cnt*th)
dt <- n_tests %>% filter(!(PATIENT_ID %in% tmp$PATIENT_ID))
fl_cols <- mean(colSums(sapply(dt[,-(1:8)], is.na)) == 0)
pnts <- length(unique(dt$PATIENT_ID)) / length(unique(n_tests$PATIENT_ID))
cls_vec <- append(cls_vec, fl_cols)
pts_vec <- append(pts_vec, pnts)
}
data.frame(threshold=th_seq, full_cols=cls_vec, patients_left=pts_vec)
}
to_plot <- veryfication_df(patient_last_results, ncol(df)-8, seq(0.5, 1.0, by=.025))
ggplot(to_plot, aes(x=threshold)) +
geom_line(aes(y=full_cols, color="Full columns")) +
geom_line(aes(y=patients_left, color="Patients left")) +
labs(color="Legend") +
geom_vline(xintercept=0.675, linetype = "longdash")
TODO
W zbiorze danych pozostawimy tylko tych, którzy mają ponad 0.675 wykonanych rodzajów badań.
threshold <- 0.675
patients_to_be_removed <- patient_last_results %>%
filter(number_of_tests < (ncol(df)-7) * threshold) %>%
select(PATIENT_ID)
df <- df %>% filter(PATIENT_ID %!in% patients_to_be_removed$PATIENT_ID)
Po usunięciu rekordów pacjentów, u których przeprowadzono mało badań, pozostało 354 unikalnych pacjentów.
knitr::kable(
df %>%
group_by(PATIENT_ID) %>%
fill(everything()) %>%
summarise_all(last) %>%
ungroup() %>%
select(-(1:7)) %>%
pivot_longer(everything()) %>%
mutate(value=is.na(value)) %>%
group_by(name) %>%
summarize(value=sum(value), .groups="drop") %>%
arrange(value)
)
| name | value |
|---|---|
| AAT | 0 |
| albumin | 0 |
| ALP | 0 |
| ALT | 0 |
| basophil_count | 0 |
| basophil_percent | 0 |
| creatinine | 0 |
| direct_bilirubin | 0 |
| eGFR | 0 |
| eosinophil_count | 0 |
| eosinophils_percent | 0 |
| gamma_GT | 0 |
| globulin | 0 |
| HCO3_ | 0 |
| hematocrit | 0 |
| hemoglobin | 0 |
| LDH | 0 |
| lymphocyte_count | 0 |
| lymphocyte_percent | 0 |
| MCH | 0 |
| MCHC | 0 |
| MCV | 0 |
| monocytes_count | 0 |
| monocytes_percent | 0 |
| neutrophils_count | 0 |
| neutrophils_percent | 0 |
| platelet_count | 0 |
| RBC_count | 0 |
| total_bilirubin | 0 |
| total_cholesterol | 0 |
| total_protein | 0 |
| urea | 0 |
| uric_acid | 0 |
| WBC_count | 0 |
| indirect_bilirubin | 1 |
| calcium | 2 |
| serum_chloride | 2 |
| serum_potassium | 2 |
| serum_sodium | 2 |
| corrected_calcium | 3 |
| hs_CRP | 3 |
| glucose | 5 |
| INR | 5 |
| prothrombin_activity | 5 |
| prothrombin_time | 5 |
| RCDW | 6 |
| RCDW_SD | 6 |
| MPV | 10 |
| P_LCR | 10 |
| PDW | 10 |
| thrombocytocrit | 10 |
| D_D_dimer | 14 |
| procalcitonin | 44 |
| aPPT | 59 |
| fibrinogen | 59 |
| thrombin_time | 59 |
| ESR | 69 |
| hs_cTnI | 69 |
| HBsAg | 82 |
| HCV_abs_quant | 82 |
| q_t_pallidum_abs | 82 |
| HIV_abs_quant | 83 |
| NT_proBNP | 90 |
| PH_value | 126 |
| interleukin_6 | 138 |
| 2019_nCov_detection | 140 |
| interleukin_1β | 140 |
| interleukin_2_receptor | 140 |
| interleukin_8 | 140 |
| TNF_alfa | 140 |
| interleukin_10 | 141 |
| ferritin | 143 |
| antithrombin | 153 |
| FDPs | 153 |
df <- df %>%
mutate(DAYS_TO_OUTCOME = as.numeric(difftime(
floor_date(DISCHARGE_TIME, "1 day"),
floor_date(RE_DATE, "1 day"), unit="days"
), unit="days"), .before=OUTCOME)
df %>% summarise(DAYS_TO_OUTCOME = quantile(DAYS_TO_OUTCOME, c(0, 0.25, 0.5, 0.75,1)), quantile = c("min", "1st", "2nd", "3rd", "max"))
## # A tibble: 5 x 2
## DAYS_TO_OUTCOME quantile
## <dbl> <chr>
## 1 -2 min
## 2 3 1st
## 3 7 2nd
## 4 13 3rd
## 5 35 max
ggplot(df, aes(DAYS_TO_OUTCOME)) + geom_histogram(binwidth = 1)
df <- df %>%
mutate(RE_DATE = floor_date(RE_DATE, '1 day')) %>%
group_by(PATIENT_ID, RE_DATE) %>%
fill(everything()) %>%
summarise_all(last) %>%
ungroup() %>%
rowwise() %>%
mutate(number_of_tests=sum(!is.na(c_across(-(1:8)))), .after=OUTCOME)
ggplot(df, aes(number_of_tests)) + geom_histogram(binwidth = 1)
df <- df %>% select(-number_of_tests)
ggplot(
df %>% pivot_longer(hs_cTnI:creatinine) %>% filter(!is.na(value)),
aes(x=OUTCOME, y=value)
) + geom_boxplot() +
scale_x_discrete(labels = abbreviate) +
facet_wrap(name ~ ., scales="free",ncol=5)
TODO
cor_mat <- cor(
x = df %>%
mutate(isMALE = as.numeric(GENDER=="MALE"),
isCURED = as.numeric(OUTCOME=="CURED"),
has_nCOV_test = as.numeric(!is.na(`2019_nCov_detection`))) %>%
select(-c(1:8), -"2019_nCov_detection"),
use="pairwise.complete.obs"
)
cor_df = data.frame(round(cor_mat,2)) %>%
rownames_to_column() %>%
pivot_longer(-rowname, names_to="colname")
cor_plot <- ggplot(cor_df, aes(colname, rowname, fill=value)) +
geom_tile() +
scale_fill_gradient2() +
theme(axis.text.x=element_text(angle = 90, hjust = 0))
ggplotly(cor_plot)
knitr::kable(
cor_df %>% filter(colname>rowname) %>% arrange(desc(abs(value))) %>% head(30)
)
| rowname | colname | value |
|---|---|---|
| MPV | P_LCR | 0.99 |
| INR | prothrombin_time | 0.99 |
| platelet_count | thrombocytocrit | 0.98 |
| HBsAg | interleukin_6 | 0.98 |
| direct_bilirubin | total_bilirubin | 0.98 |
| lymphocyte_percent | neutrophils_percent | -0.98 |
| procalcitonin | q_t_pallidum_abs | 0.95 |
| D_D_dimer | FDPs | 0.95 |
| P_LCR | PDW | 0.95 |
| MPV | PDW | 0.94 |
| ALT | interleukin_1β | 0.93 |
| lymphocyte_count | monocytes_count | 0.88 |
| serum_chloride | serum_sodium | 0.86 |
| eosinophil_count | eosinophils_percent | 0.86 |
| AAT | interleukin_1β | 0.86 |
| hematocrit | hemoglobin | 0.84 |
| MCH | MCV | 0.82 |
| indirect_bilirubin | total_bilirubin | 0.80 |
| interleukin_6 | interleukin_8 | 0.80 |
| RCDW | RCDW_SD | 0.79 |
| aPPT | prothrombin_time | 0.79 |
| monocytes_percent | neutrophils_percent | -0.77 |
| eGFR | urea | -0.77 |
| creatinine | urea | 0.75 |
| albumin | calcium | 0.74 |
| isCURED | neutrophils_percent | -0.73 |
| isCURED | lymphocyte_percent | 0.72 |
| neutrophils_count | neutrophils_percent | 0.71 |
| albumin | total_protein | 0.70 |
| lymphocyte_percent | neutrophils_count | -0.70 |
knitr::kable(
tidy(cor.test(df$LDH, as.numeric(df$OUTCOME=="CURED")))
)
| estimate | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|
| -0.6270327 | -24.33491 | 0 | 914 | -0.6648056 | -0.5860615 | Pearson’s product-moment correlation | two.sided |
correlations <- df %>% mutate(isCURED=as.numeric(OUTCOME=="CURED")) %>% select(-c(1:2,4:8)) %>%
pivot_longer(!isCURED, names_to="attribute", values_to="x_val") %>%
filter(!is.na(x_val)) %>%
nest(data = c(x_val, isCURED)) %>%
mutate(cor_test = map(data, ~tidy(cor.test(.x$x_val, .x$isCURED)))) %>%
unnest(cor_test) %>%
select(attribute, estimate, p.value, conf.low, conf.high) %>%
arrange(desc(abs(estimate))) %>% head(15)
## Warning: Problem with `mutate()` input `cor_test`.
## ℹ the standard deviation is zero
## ℹ Input `cor_test` is `map(data, ~tidy(cor.test(.x$x_val, .x$isCURED)))`.
## Warning in cor(x, y): the standard deviation is zero
knitr::kable(
correlations
)
| attribute | estimate | p.value | conf.low | conf.high |
|---|---|---|---|---|
| neutrophils_percent | -0.7301461 | 0 | -0.7586392 | -0.6988654 |
| lymphocyte_percent | 0.7234883 | 0 | 0.6915912 | 0.7525691 |
| albumin | 0.6880859 | 0 | 0.6524294 | 0.7207028 |
| prothrombin_activity | 0.6547754 | 0 | 0.6084767 | 0.6966318 |
| hs_CRP | -0.6509896 | 0 | -0.6910468 | -0.6069460 |
| D_D_dimer | -0.6376192 | 0 | -0.6819778 | -0.5885867 |
| LDH | -0.6270327 | 0 | -0.6648056 | -0.5860615 |
| neutrophils_count | -0.6083675 | 0 | -0.6470960 | -0.5665074 |
| FDPs | -0.6042409 | 0 | -0.6689583 | -0.5304310 |
| AGE | -0.5761382 | 0 | -0.6071595 | -0.5433633 |
| calcium | 0.5703332 | 0 | 0.5257539 | 0.6117881 |
| platelet_count | 0.5419441 | 0 | 0.4952125 | 0.5855487 |
| monocytes_percent | 0.5415337 | 0 | 0.4947996 | 0.5851444 |
| eGFR | 0.4909936 | 0 | 0.4402769 | 0.5385870 |
| urea | -0.4841428 | 0 | -0.5321758 | -0.4330031 |
tseries <- df %>%
select(1:8) %>%
group_by(PATIENT_ID) %>%
summarize(across(everything(), last), .groups="drop") %>%
select(-RE_DATE, -AGE, -GENDER, -DAYS_TO_OUTCOME, -PATIENT_ID) %>%
pivot_longer(ADMISSION_TIME:DISCHARGE_TIME, names_to="TYPE", values_to="DATE") %>%
mutate(DATE = floor_date(DATE, "days")) %>%
mutate(TYPE = str_replace(TYPE, "_TIME", "")) %>%
count(TYPE, DATE, OUTCOME) %>%
rename("COUNT" = "n")
tplot <- ggplot(tseries, aes(x=DATE, y=COUNT, fill=OUTCOME)) +
geom_col() +
scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b %Y") +
facet_grid(TYPE ~ .)
ggplotly(tplot)
unique(c(.packages(), loadedNamespaces()))
## [1] "plotly" "tibble" "broom" "purrr" "ggplot2"
## [6] "stringr" "lubridate" "dplyr" "tidyr" "arules"
## [11] "Matrix" "httr" "readxl" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
## [21] "Rcpp" "highr" "cellranger" "compiler" "pillar"
## [26] "tools" "digest" "viridisLite" "jsonlite" "evaluate"
## [31] "lifecycle" "gtable" "lattice" "pkgconfig" "rlang"
## [36] "cli" "crosstalk" "curl" "yaml" "xfun"
## [41] "withr" "knitr" "htmlwidgets" "generics" "vctrs"
## [46] "grid" "tidyselect" "data.table" "glue" "R6"
## [51] "fansi" "rmarkdown" "farver" "magrittr" "backports"
## [56] "scales" "htmltools" "ellipsis" "assertthat" "colorspace"
## [61] "labeling" "utf8" "stringi" "lazyeval" "munsell"
## [66] "crayon"